arXiv:1503.07560vl [math.NA] 25 Mar 2015 


A REDUCED BASIS METHOD FOR THE 
HAMILTON-JACOBI-BELLMAN EQUATION WITH APPLICATION 
TO THE EUROPEAN UNION EMISSION TRADING SCHEME* 

SEBASTIAN STECK^ AND KARSTEN URBAN§ 

Abstract. This paper draws on two sources of motivation: (1) The European Union Emission 
Trading Scheme (EU-ETS) aims at limiting the overall emissions of greenhouse gases. The optimal 
abatement strategy of companies for the use of emission permits can be described as the viscosity 
solution of a Hamilton-Jacobi-Bellman (HJB) equation. It is a question of general interest, how 
regulatory constraints can be set within the EU-ETS in order to reach certain political goals such 
as a good balance of emission reduction and economical growth. Such regulatory constraints can be 
modeled as parameters within the HJB equation. 

(2) The EU-ETS is just one example where one is interested in solving a parameterized HJB 
equation often for different values of the parameters (e.g. to optimize their values with respect to a 
given target functional). The Reduced Basis Method (RBM) is by now a well-established numerical 
method to efficiently solve parameterized partial differential equations. However, to the best of our 
knowledge, an RBM for the HJB equation is not known so far and of (mathematical) interest by its 
own, since the HJB equation is of hyperbolic type which is in general a nontrivial task for model 
reduction. 

We analyze and realize a RBM for the HJB equation. In particular, we construct an online- 
efficient error estimator for this nonlinear problem using the Brezzi-Rapaz-Raviart (RBB) theory as 
well as numerical algorithms for the involved parameter-dependent constants. Numerical experiments 
are presented. 

Key words. Reduced Basis Method, Hamilton-Jacobi-Bellmann equation, emission trading 
system 
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1. Introduction. A driving source of motivation for the investigations reported 
in this paper is the European Union Emission Trading Scheme (EU-ETS) that has 
been invented in order to limit the emission of greenhouse gases. Within this EU- 
ETS, a limited amount of emission permits is issued and each pollutant needs to cover 
its emissions with sufficient permits. Both from the environmental and ecological as 
well as from the economical point of view it is important to control the EU-ETS in 
such a way that certain desired political effects are reached. As a simple example, the 
number of permits should limit the emission of global warming gases without leading 
to a severe economical crises. Thus, we are interested in investigating the effect of 
several different regulatory constraints (i.e., a parameter study from a mathematical 
point of view) as well as trying to find regulatory strategies in order to reach certain 
goals (i.e., realtime optimal control). 

Erom a mathematical point of view, this means that the same model has to be 
solved for a variety of parameters, here different regulatory constraints. We describe 
these constraints in terms of parameters /r G X>, where T> C K.^ is the set of all 
possible parameter values. Moreover, a full mathematical model of the EU-ETS is in 
general complex so that the numerical simulation is costly and parameter studies as 
well as realtime optimal control is not feasible. Hence, we suggest to use the Reduced 
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Basis Method (RBM), a model reduction technique that uses a possible costly offline 
phase in order to computationally construct a reduced system, which is then used in 
the multi-query (parameter study) or realtime (optimal control) context in order to 
produce numerical approximations for various parameter values highly efficient and 
with mathematical certification in terms of a posteriori error control. 

It turns out that an optimal abatement strategy concerning the use of emission 
permits can be described by the Hamilton-Jacobi-Bellman (HJB) equation, parame¬ 
terized by the regulatory constraints. This brings us to the second source of motivation 
for this paper, which is of mathematical nature and independent of the specific ap¬ 
plication of the EU-ETS. In fact, to the best of our knowledge, an RBM for the HJB 
equation is not known so far and of interest by its own since the HJB equation is of 
hyperbolic type which is in general a nontrivial task for model reduction. 

It is the aim of this paper to construct, analyze and realize a RBM for the HJB 
equation with the specific application of the EU-ETS. The remainder of this paper 
is organized as follows. In Section 2, we collect the required preliminaries on the 
mathematical model for the EU-ETS yielding the HJB equation. Section 3 is devoted 
to the discretization including an error analysis. Since we are facing a nonlinear 
problem (in the specific case of emission trading a quadratic problem, see also Remark 
2.1 below), we use the Brezzi-Rappaz-Raviart (BRR) theory in order to ensure well- 
posedness and error control. The RBM for the HJB equation is introduced in Section 
4 and finally, in Section 5 we present results of some numerical experiments. 

2. Preliminaries. In this section, we collect some preliminaries. 

2.1. Mathematical Model for the European Union Emission Trading 
System (EU-ETS). We start by introducing a mathematical model for the Euro¬ 
pean Union Emission Trading System (EU-ETS) and show that the market equilib¬ 
rium can be described in terms of a Hamilton-Jacobi-Bellman (HJB) equation. First, 
the emission trading system is organized in trading periods, but for simplicity we may 
reduce ourselves to one period only. It can be shown under reasonable assumptions, 
that the market equilibrium for one single trading period [0,T] can be characterized 
by the fact that the sum of the costs of all market participants is minimal, [8]. 

Let Yr be a stochastic process, r S [0, T], which describes the amount of uncovered 
emissions, sometimes also called state. For a set of d companies whose greenhouse 
gas emissions are considered, the values of U- are thus taken in The valued 
stochastic control tTt describes the additional abatement of emissions compared to 
the so-called business as usual strategy. Hence, an optimal abatement strategy tt = 
(7rT-)^g[Q should minimize the expected abatement costs, i.e., the cost functional 



( 2 . 1 ) 


where f'^ denotes the running abatement cost using strategy tt and the function h 


models the penalty to be paid at the end of the trading period. 


Remark 2.1. For later reference, we note that in the specific case of the EU-ETS, 
the dependency of f'^ with respect to tt is quadratic. 

A standard stochastic model for the amount of uncovered emissions U- reads 



where WV is a d-dimensional Wiener process and b^, are coefficients for drift and 
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volatility, respectively, such that is linear w.r.t. the control tt as well as 
is linear in tt. 

The transition to a partial differential equation (PDE) in terms of a HJB equation 
is then done by introducing new variables (t, x) for time and state and to modify ( 2 . 2 ) 
to 


(2.3) 


(IYt = b'^{T,Yr)dT + {T,Yr)dWT, 


T e {t,T], 


Yt = X, 


i.e., the initial time t and the initial state x at t = t are the new variables. Moreover, 
the stochastic control tt is replaced by a determinstic (but state-dependend) function 
7 : {t,x) I—>• and accordingly, the cost functional in ( 2 . 1 ) reads 


(2.4) 


J{t, x\ 7 ) := IE 


P {TPr)dr + h{YT) 


where the dependency of x is implicit via Yt = YT{t,x) by (2.3). In order to derive a 
strategy that yields minimal cost, one needs to solve a stochastic optimization problem 
whose solution is the value function, i.e., for x € 

(2.5) u{t,x) = ini J{t,x]p Vt€[0,T), u{T,x) = h{x), 

7er 


where E C Too([0, T] x Mf") is a suitable set of admissible controls. It is well-known 
that the value function is a viscosity solution of the Hamilton-Jacobi-Bellman (HJB) 
equation (see e.g. [17, Theorem IV.5.2]) 

(2.6) dtu{t, x) + sup I pT:{a'^{a'^p V^u(t, x)) + ■ Vu{t,x) — P{t, a:)| = 0, 

7er 12 J 

for all (t, x) £ [0, T] x B, where denotes the Hessian of u. We consider a bounded 
domain B C for the state. The reason is twofold: (1) The limit of the price 
for emissions exceeding the available permits is the penalty set by the regulating 
authorities; (2) If, on the other hand, there is a vast excess of permits, no emissions 
will be saved and the value of the permits is only determined by their terminal value 
at T. Of course, using a bounded domain 11 C R"^ implies the need to set appropriate 
boundary conditions. Since the price for the permits is determined by the derivative of 
the value function u w.r.t. the need for permits x, a corresponding Neumann boundary 
conditions are appropriate, see ( 2 . 8 b) below. 

Last, but not least, we model the appearance of parameters n £ 1) C R'^ as 
already introduced in Section 1 . This means that basically all quantities may be 
/^-dependent, e.g. P{p), b'^{p, J{fj,;t,x;p and the value function u{fi) = 

up; t, x). 

2.2. Hamilton-Jacobi-Bellman (HJB) Equation. Let T > 0 be some final 
time, I := [0,T] the time interval, B C R'^ be a bounded domain and denote by 
VIt := I X n the time-space domain. The parameters are denoted by /x G U, where 
V C R^ is the parameter space. 

Next, let r C Loo(f2T;R'^) be a separable, complete metric space, the space of 
admissible controls, such that the mapping 

r ^ c{n) X c{n, r'^) x cp) x c'(b), r 9 7 {a^p),b^p),c'p),p{P)) 
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is continuous for any parameter fj, G T>. We define the linear parameterized PDE 
(PPDE) operator in space by 

u I—>■ u) := —a^(/i) Ait + b~*(fj-) ■ Vit + fj, GV. 

In the application of the EU-ETS described in §2.1 above, this operator takes the 
form := V'^u) + 6 ^(/i) • Vm. Note, that A'^{^]u) is linear 

in It, whereas the dependency on the control 7 is in general not linear. Einally, we 
define the nonlinear parameterized Hamilton-type operator by 

(2.7) u) := sup{A'f{^^■, u) - 

7er 

as well as the linear space-time differential operator it) := dtu + A^{^\ it), which 

is linear in it, but not self-adjoint. The range of this operator should also contain 
boundary and terminal conditions included in the following set of Hamilton-Jacobi- 
Bellman (HJB) equations (compare (2.6)) 

( 2 . 8 a) dtu+ 'H{y.\u) = 0, in Hy, 

d 

(2.8b) —u = i/>, on = (0, T) X 9H, 

( 2 . 8 c) w(T) = ut, on H, 

where '0 is a suitable function modeling the Neumann truncation boundary conditions 
and Ut G H^(fl) denotes the terminal condition (which, in particular needs to be 
compatible with 0 on {T} x dfl). For the correct interpretation of the subsequent 
discussion, it is worthwhile to detail (2.8a) in combination with (2.7): 

(2.9) dtu{t,x)+ sup {A^(*’"')(^;it( 0 a:)) -/'^(*’'')(^)} = 0, x) G 

'y{t,x)GR^ 

This shows that the control space is P := Loo(f^T; R'^), or an appropriate subspace 
of the latter. We can also write (2.8) in the following form 

( 2 . 10 ) sup{L'’'(^; 11 ) — g'’'(^)} = 0 on Ht- 

Given a value function, i.e., a solution u* G U (where U is an appropriate solution 
space, e.g. H^{I; iJ“^(H))nL 2 (/; i7Q(H))) of ( 2 . 8 ), the corresponding optimal control 
7 * S P is given by 

7 * = argsup{L'^(Ai;u*)-g^(/r)}. 

7Gr 

The pair of optimal value function and optimal control is also written as x* = 
gX:=T xU. 

Well-posedness. It is well-known that well-posedness of the HJB-equation (2.8) 
is ensured under reasonable assumptions. In fact, if A'^, /'>' and ut are uniformly 
continuous and uniformly Lipschitz continuous with respect to the state x G as 
well as bounded at the state a; = 0, it was proven e.g. in [17, Theorem IV.6.1] that 
the HJB equation (2.8) admits at most one viscosity solution such that there exists 
constant K gM. with 

(2.11a) |M(t,a:)| 0 i4r(l-I-|a:|) 

( 2 . 11 b) \u(t,x) — u{t,x)\ 0 i4r(|a; — a;| -I- (1 -I- min{|a;|, |a:|})|t — 


4 


for all {t,x), {i,x) G fir- 

Furthermore, the value function u in (2.5) is a viscosity solution of the HJB 
equation (2.8), which satisfies (2.11) [17, Theorem IV.5.2 & Proposition IV.3.1]. Thus, 
the value function (2.5) is the unique solution of the HJB equation. 

Remark 2.2. It should be noted that we do not have a linear-quadratic problem 
even though /'>' is quadratic in 7 , the HJB remains nonlinear and we cannot expect 
the availability of a solution formula or even a smooth solution, [17]. 

3. Discretization. We now describe the essentials of a numerical discretization. 
We start by a possibly high-dimensional model that is assumed to reflect the true 
model sufficiently well, thus called ‘truth’ discretization. This will later be the basis 
for model reduction. 

3.1. ‘Truth’ discretization. Recall from (2.9) that the supremum is taken 
pointwise. This also motivates that most discretizations are pointwise. In order to 
describe such schemes, let Zi := {ti,Xi) G By, i = 1,... ,JV, be a set of points in the 
time-space domain, where W ^ 1 is assumed to be ‘large’, in particular large enough 
to represent the main characteristics of the continuous problem ( 2 . 10 ) as well as ‘too 
large’ for multi-query or realtime simulations. Hence, we are looking for a discrete 
approximation = (u*’^(/r))i=i,,,,,AA S of the value function u* G U at 

the points Zi = (ti,Xi). This means that a pointwise discretization of (2.10) amounts 
solving an optimization problem for each point, i.e., we need to determine a component 
of an approximation to 7 * G T for each i. We denote such an approximation of 
'y*(ti,Xi) by 7 *’^(^). If we abbreviate the pointwise evaluation of the operator and 
right-hand side, respectively, as 

•) : ^ K, ■) := 

/^■'^•(m) G K, /^’^“(m) := [r(*‘’"‘)(M)](i.,^.), 

we obtain the discretized optimization problem of dimension A/”: 

(3.1) Find M^(^) G : max{L-^’T'*(^;u^(^)) -/•^’T'*(^)} = 0. VI < i < W. 

7iGR'^ 

The corresponding solution is denoted by u*'^[p) = (u*’^(/r))i=i_..._// G which 
is a discrete approximation of the value function u* G U. For the optimal control, we 
set (with the solution u*'-^{p) of (3.1)) 7 *’^(Ai) = {p))i=i,...,M G defined 

for each 1 < i < A/" by 

(3.2) 9 ■■= arg {p)) - f^'^'{p)}, 1 < i < Af. 

7ieR‘^ 


Remark 3.1. The above described model yields F = Loo(f2T; R”*) or an appropri¬ 
ate subset. Such a subset would occur if control constraints would appear. In that case, 
R^^ would be replaced by some E C R'^. Correspondingly, R"^-^ in the discretization 
would have to be replaced by E^. We emphasize that all subsequent findings remain 
true in this case. 

In order to simply notation, we collect all single optimization problems into one 
system by setting for 7-^ G G R'^ 
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which means that ■) : x R)-^ —>■ R-^ and /^’'(/^) : R'^'^ —>■ R"^. Hence, for 

the pair 


(7*’^(/r),M*’^(/r)) G := (R^ x R)^ = r('^+i)^ 

of optimal (discrete) control and optimal (discrete) value function, we get that 

L^’T'*''^('^)(Ai;u*’^(/r)) - /^’^*’^('^)(Ai) = 0 in R^. 

3.2. A nonlinear system for control and state. If the function 
L^’'>'^(/r,u*’^(/r)) — ip) is in C'^(R‘^;R)^, we can consider the Frechet deriva¬ 

tive w.r.t. the control variable 7 at some <5^ G R*^^ 

G L(R^,R)^, 9T,[/^’'*^(Ai)] G L(R^R)^.i 


Recall that in the case of the EU-ETS, p is a quadratic function of the control 7 so 
that the assumed differentiability in fact holds. Then, the optimal control 7*’^(/r) G 
R"^^ is a critical point of this mapping, i.e. 

(3.3) a^[L^’'>'*’^('^)(Ai;u*’^(Ai))-/^’^*’^('^HM)](<5^) = 0, V<5^ G R'^^, 

which is a nonlinear problem for (even though (3.3) is linear in the ‘test 

function’ 5^). We define the composite function [p : X-^ := (R'^ x R)-^ —>■ 
L(R‘^,R)^ X R^ =: as {x^ = ,u^) G X^) 


(3.4) g^iPix^) := 


P;u^) - {P]\ ^ 


In this notation, the ‘truth’ control/state-solution x^'*{p := p*{p, u*(p) G 
X-^ is characterized by 


(3.5) 


g^iPP^’PP) = g^{PP*’^{P,u*’^P)) = 0 in Y^, 


i.e., this equation is to be understood in L(R‘‘*,R)^ x R-^ = Y-^. 

For later purpose, we determine the Frechet derivatives of fy^(/r), in case of 
their existence (which is obviously guaranteed for the EU-ETS case), of course. 
Let = p^,u^) G X^. Then, we have D{g^{P){x^) G L(X^,Y^), i.e., 
{D{g^{p)p^)){x^) G Y^ for x^ = pp'p-^) G X-^ and get 


{D{g^{P){x^)){i^) 


(dpgp{p{xP))pp + dpgp{p{xP)){up\ 
\dpgp{ppp)pp + dpgpp)pp){up) 


f d, p, [l^Ai^-. (m)] ) (7^) + a„ {d^ [l^Ap\ P)]) pp] 

[ dpL^^^^P;up-A^^P))pp+dpL^’^^P;up-A^^iP)iup ) 


[ d, if,; up - iP] PA + (/r; u^) ) ' 


^We denote by L{X, Y) the space of continuous, linear mappings from a normed space X to a 
normed space Y. 
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Remark 3.2. One possible numerical method to determine a solution of (3.1) 
is the so-called policy iteration algorithm, also called Howard’s algorithm, [3], which 
reads as follows for an initial value For k = 0,1, 2,... do 

(3.7a) = arg max {L^’^(/x; «('=)) - /^’^(a^)}, 

(3.7b) findu^’^+^^ (^i; (m)- 

Under appropriate conditions, this algorithm converges and for the limits, we have 
—>■ u*’^{pl) (the solution of (3.1),) —>• 7 *’^(/r) as fc —>• oo (defined in (3.2)/ 

For later use in deriving error estimates, we collect conditions that ensure Lip- 
schitz continuity of the Frechet derivative DQ-^{p,). 

Lemma 3.3. If the estimates 

(3-8a) < gg IIt/^ — T^IIgdAf, 

(3.8b) - L^’'^^w^)||i(KdAr^RW) < g^Wxf - ||xw, 

(3.8c) - L-^’^^w^)||L(Rdw^L(Rdw_KAf)) < 02 \\xf “ IIx-a, 

hold for constants < oo, k = 0 , 1 , 2, Xi = (yi, Ui) G X-^, i = 1 , 2 and 

(3.9a) ||5^(/^’'^^ -/^■^^)||l(«^a.,r..) < g{\hf 

(3.9b) \\9^{f'^''^^ ~ )llL(R‘i-^,L(R<iJv^RAr)) < ~ 

for constants gj, < oo, k = 1,2, then D{Q-'^{ij)) is Lipschitz continuous, i.e., 

(3.10) \\D{g^ip))ixf) - 7^(l?^(A^))(xf )|U(X^,YW) < g\\xf - x^||x 

with g = 00 + 2 gf 02 + g( + g( ■ 

Proof The proof is more or less standard using the representation (3.6) of 
D{Q^{pi)) and triangle inequalities several times. It only remains to note that using 
= u-^ = U 2 in (3.8b) implies 'a'^)llL(RdAf^RAf) < giW'^'i — 

7 ^ 11 xJX, which results in the estimate || 57 (L^’'’'^-L^’'’'^)||L(RXf,L(R<i-^,KAf)) < g{\\'^^- 
7 ^||xw, which is needed in the estimate (3.10). □ 

Remark 3.4. Recall that all assumptions in Lemma 3.3 are satisfied for the 
EU-ETS. In fact, the operator L'^ is linear in 7 and is a quadratic function of j. 

3.3. Well-posedness and error control. We now present some well-known 
ingredients of the theory developed by Brezzi, Rappaz and Raviart, known as BRR 
theory, [4, 5, 6 ], which provides us with results concerning existence and uniqueness 
as well as error control for the nonlinear system (3.5). In order to streamline the 
notation, we consider a general generic framework. To this end, let X, Y be two 
finite-dimensional spaces normed by || ■ ||x: and || • ||v, respectively. We consider a 
nonlinear mapping G : X ^ Y and seek for a solution x* G X of the problem G{x) = 0 
in Y (i.e., a general form of (3.5)). 

Next, we assume that the Frechet derivative DG{x) G L{X, Y) exists for all x G X 
and that the inverse also exists. Then, we define for some fixed x € X the mapping 

IIs(x) :=x- {DG{x))-^Gix), 


(3.11) 


H^:X^X, 
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which is obviously of quasi-Newton type. Finally, we set 

(3.12) /?£:= ||pG(x))-i||-(V,x) 
and assume that 0 < /3g < oo. 

Lemma 3.5. Let the Frechet derivative DG{x) € L{X,Y) exist for all x S X, he 
invertible and Lipschitz continuous, i.e., there exists a constant g > 0 such that 

(3.13) \\DG{xi) - DG{x 2 )\\l{x,y) < 9\\xi - X 2 \\x 
for all Xi^X 2 € X Then, for any Xi,X 2 € X, we have 

pi 

\\Hg,{xi) - Hs;{x 2 )\\x < -f \\xi - X 2 \\x / \\x - {X 2 + t{xi - a; 2 ))|U dt. 

Px Jo 


Proof. The proof follows standard lines starting with the fundamental theorem 
of calculus G{xi) = G{x 2 ) + Jq DG{x 2 + t[xi — X 2 )){xi — X 2 ) dt. Then, we get 

Hs:{xi) - Hs{x 2 ) =Xi-X 2 - {DG{x)) ^ {G{xi) - G{x 2 )) 

= {DG{x)y\{DG{x)){xi - X 2 ) - {G{xi) - G(a; 2 ))} 

= [DG{x)) ^^[DG{x)){xi - X 2 ) - J DG[x 2 + t{xi - X 2 )){xi - X 2 ) df 

pi 

= {DG{x)) W {DG{x) - DG{x 2 + t{xi - X 2 ))}ixi - X 2 ) dt. 

Jo 

Next, we use standard estimates to obtain 


1 r 

\\Hs,{xi) - Hs{x2)\\x < ^ / \\{DG{x) - DG[x2+t{xi - X2))}{xi - X2)\\Ydi 

Px Jo 

Q /* ^ 

<-f\\xi-X 2 \\x / \\X - {X2 +t{xi - X2))\\x dt, 

Px Jo 

by (3.13), which proves the claim. □ 

Lemma 3.6. Let the assumptions of Lemma 3.5 hold. Then, the mapping is 
a contraction on B~^{x) ■.= {x & X ■. \\x — x\\x < 7 } if 1 < 7contr. := 

Proof. Let xi,X2 S B.y{x). Then, we have \\x — (x2 -I- t{xi — 2:2))||a < 7 and the 
assertion follows immediately by Lemma 3.5. □ 

Lemma 3.7. Let the assumptions of Lemma 3.5 hold. Then, the mapping Tlx '■ 
Bj{x) —>■ Bj{x) is a self-map for all x G X with 

(3.14) r{x):=^JGix)\\Y<l, 

and 


Px 


7 e [7min,7max] := — 1 - - t{x), 1 + ^/T 

Q L 



(3.15) 






Proof. We start by the simple identity H^{x) — x = Hx{x) — H^{x) + H^{x)—x = 
Hx{x) — Hx{x) — {DG{x))~^G{x), which holds for for any x € X. Then, 

\\Hx{x) - x\\x < \\Hxix) - Hxix)\\x + /3j ^||G(x)||y. 


If a: G B~f(x), we can further estimate the first term by Lemma 3.5 and obtain 


\\Hs{x) - 


x\\x < 


El 
Px Jo 


t\\x - x\\x dt + I3 ^'^\\G{x)\\y < 1^+/3s^l|G(x)|h 


This latter term is less than 7 if and only if 7 ^ — ^^7 + | ||G(x) ||y < 0, which in turn 
is valid for 7 e [ 7 min, 7 max] and ||G(x)||y < (2p)“i/3?. □ 

Summarizing the above findings, we get the following result. 

Proposition 3.8. Let the Frechet derivative DG{x) € L{X,Y) exist for all 
X € X, be invertible and Lipschitz continuous with constant g. Let x € X be given 
such that (3.14) holds. Then, there exists a unique fixed-point x* S Bj{x) of Hx for 

all 7 € [7min j 7 contr.) ■ 

Proof. The proof follows from Banach fixed-point theorem in view of Lemma 3.5 
and 3.7. □ 


Proposition 3.8 yields a well-posedness result, but at the same time also provides 
us with an error estimate as we shall see next. 

Corollary 3.9. Under the assumptions of Proposition 3.8, the estimate 


(3.16) Ik* - x\\x < —(1 - \/l - r(a;)) 

Q 

holds for X G X satisfying (3.14). 

Proof. The mapping Hx has a unique fixed-point x* in B-f[x) for 7 = 7min. □ 

Remark 3.10. The following observations are potentially crucial for the numer¬ 
ical realization: 

(a) Note, that the quantity t{x) is an a posterori indicator provided that g (or an 
estimate) is available and /3x (or an estimate) is computable. In fact, ||G(a;)||y 
is the computable residual of the nonlinear equation and thus (3.14) can be 
verified a posteriori. We will later use this observation to obtain an error 
estimate for a numerical approximation x of the solution of the nonlinear 
problem G{x) = 0. 

(b) Obviously, the assumptions on the Frechet derivative DG only need to hold 
in a neighborhood of x. 

4. The Reduced Basis Method (RBM). We now introduce the Reduced 
Basis Method (RBM) for the numerical solution of the parameterized HJB equation 
and start by reviewing the standard RB setting. The main idea is to select in an 
offline phase so-called snapshot parameters 


Sn '■= {mi, • ■ ■, Ma}, 

and compute the corresponding snapshots Xi := x*'-^{pLi) = {pLi)) G 

X-^ as the solution of (3.1), (3.2), or in other terms, (3.5) in Y-^. The precise selection 
of the snapshots will be explained later. Then, we define the RB space as Xx '.= 
spanja:)^ : f = 1,..., N}, assuming that N <^Af. 
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An RB approximation x*j^{fi) € of G X-^, € T> \ Sn, is then 

computed by the solution of 


g{fi){xl,{fi)) =0 inYNin), 

where C is some test space which is possibly ^-dependent and such that 

the reduced problem is stable. 

The aim is to use the RBB-theory for developing an a posteriori error estimate 
for the RB-approximation. To this end, we fix some x € X-^ (to be determined below) 
and define for any = ( 7 ^, € X-^ the mapping 

(4.1) H^{fJ,)ix^) := - ((Dg^(/r))(x^))-ig^(/x)(x^), H^ip) : X^ ^ X^. 


The idea is to use x = x%{^) = ( 7 ^(^), u^(/r)) and set analogously to (3.12) 
(4.2) PisfitJ-) :=/3^j,(^)(Ai) := \\{{Dg^{ti)){x*ff{^)))-^\\ 

L{Y^,X^) 

\\iDg^it,){xU^mx^)\W^r 


We mention [1, 2] for a POD-based model reduction approach of the HJB-equation, 
where the reduction is performed with respect to time. 

4.1. Computation of a lower inf-sup bound. To obtain a lower bound for 
the inf-sup constant /3Ar(/r), we follow an idea presented in [16]. We start by detailing 
the computation of a lower bound for /3Ar(/r). Since = L(R'^,K)^ x is a 
Hilbert space, we get 


/3Af(M) 




For any given x,x € X^, we define the supremizer s{p] x, x) € Y^ as 

((Dg■^(p)(x))(x),y-^)Y^f 


(4.3) 
so that 


s{fj,; X, x) := axg sup 


Y^r 


^Niy) = inf 


> inf 


{{Dg^{fi){x*jY{p.))){x^), s{fi; xllip.), x^))yJY 
{{ng^ {p){x*j^{p))){x^), s{^i] x^))y^j- 


where s{p]x^) G Y-^ is arbitrary and will be chosen later. Then, we obtain 


Pn{p-) > 


||s(/i;x^)||YV 


inf 


jjx^l 

(4.4) =: 


inf 

gx^ 


{{Dg^{^l){xN{fJ.)))ix^),sifJ.■,x^))Y^r 




where ,3£®'"®(/r), /3 lb'“(m) nre lower bounds to be computed offline and online, re¬ 
spectively. 

We determine offline a small set of i? G N so called anchor parameters Sr := 
{fLi,..., Pr} C X>. Online, for a given /r G H, we determine the ‘closest’ anchor 
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parameter S Sji defined by minimizing — p,\ over /2 G Sn. Since we choose 
s{fj,;x^) := it holds /3 lb™'^(a^) = PNiP-ifJ-))- Thus, in the offline- 

phase r = 1,..., i?, are precomputed solving generalized eigenvalue problems. 

To determine the anchor points p,r, we use a Greedy-type method detailed in Algo¬ 
rithm 1 based upon a training set C T). Of course, the constant i in line 2 

of Algorithm 1 could be verified. This choice is motivated by (4.5) since Algorithm 1 
ensures that min /3 ?r™®(u) > ^ whose relevance will be described next. 

^^^train 


Algorithm 1 Greedy selection of anchor points. 

1 : choose fii €V arbitrarily, N 1, Sr := {fii}, compute /3 lb’™(Mi) 
2: while min /3 lb™'’(m) ^ -5 do 

..^Tranchor ^ 

^^“train 

3: p.N +1 ^ argmin ,d£B“®(/r) 

^'^^train 

4: Sr ^ SrU {p-N+l}, Compute /3£®‘“(/iAr+i) 

5: N^N+1 

6 ; end while 


The online part /?lb™'^(/^) bound relies on a separation of DQ-^{p) with 

respect to the parameter and is computed by the well-known Successive Gonstraint 
Method (SGM) [12], which will be discussed below. 

Altogether, we obtain an inf-sup lower bound 

(4.5) /3LB(^) := /3£®“(/2(^)) “(m) > i/3£f “(m)- 

In order to compute an upper bound for the indicator in (3.14), we define 

(4.6) rjv(Ai) := ll^^(M)(a:Ar(M))llYAA 

and the corresponding upper bound 

(4.7) r™( a^) := \\g^{ti){x*M)\\Y^■ 

Using (3.16), the error can thus be estimated with 

(4.8) A„(m) A-^riO). 

4.2. Offline/online-separation. For an efficient separation of the required cal¬ 
culations into a possibly costly offline and a highly efficient online phase, one usually 
requires a separation of the parameter and other types of variables, sometimes also 
called affine decomposition. Here, this means 


Ql 

(4.9) L'^(/i;M) = ^d^(/x)L^(u), 

9=1 


Qf 

r(/r)=^d7(/r)/;. 


9=1 
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with functions : X> —>■ R. This separability in the parameter also transfers to 

their discretized variants and to DQ so that (3.6) reads as 




9=1 


Qf 


9=1 






(4.10) 


Ql Qf 

r ff, _ AC. „ fV 


E< (/i)l^e,^(xE(i^)-E< (/^)^^9^(i^)- 

g=l 9=1 


By inserting (4.10) and the representation XN^fj,) = in the definition 

(4.4) of /3 lb‘™(/i) we obtain: 


/3“(m)= inf 

a;AC gx^c 


{(Ht/^(^)(a;jv(M)))(ai^), s{p; x^))yac 


l|5(/i;a:-^)|lYV 


Qi Af 


(4.11a) 

(4.11b) 


+ 




n=l 

Qf 

-E 

9=1 


{DGq i^„)ixr ), s{p] X^ ))Yf^ 
I|s(m;2;^)IIyac 


(ll) 

(n)\ 

{Dgf{x^),s{g.]x-^))Yff 

1 VPJ 

^q VMJ J 

l|s(/i;a;^)llYAr 


In order to apply the Successive Constraint Method, we summarize the terms in 
(4.11a) and (4.11b) in T(^;z) 

Ql n 


:=E E {n)x’^{n)) z„+(q_i)jv 

q—1 n—1 

Qf 

+ E ZQLN+q- 


9=1 


Since we have chosen s(/x;x^) = (DQ-^{p){xNip))){x^), we get the representation 
/3 £b““(ai) = 1 + inf r(/i;z) with 




Y _ / j, p mQi,lV+Qj- . q p _ {^Sqi^n){x^),s{fJ,;X^))Yf^ 


^QLN+q — 


{DG l{x^),s{^i\x^))Yff 
l|s(/i;a:^)llYAf 


To obtain a lower bound for inf^g Zn z) we use the Successive Constraint Method 

from [12] and use an appropriate superset of Z^. 

Next, we need the Lipschitz-constant g in (3.10), which in the general case could 
also be /i-dependent. In the simpler case g ^ g{p), this constant can be computed up 
to numerical precision by solving a generalized eigenvalue problem offline. 
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Finally, it remains to compute the residual ||0^(/r)(a;^(/i))||y, which depends 
on the problem at hand and also requires the above mentioned separation. We will 
detail this for a specific numerical example in Section 4 below. Combining all this, 
the inf-sup lower bound in (4.5), the indicator in (4.6) and the error 

estimator can be computed online-efiicient. 

4.3. Greedy Algorithm. Now, we describe the computation of the snapshot 
parameters Sn. To this end, we use a standard Greedy method over the error es¬ 
timator An{p) based upon a training set Strain C T). The scheme is displayed in 
Algorithm 2. As we have seen, the error estimator AAr(/j,) is only meaningful if the 
indicator TAr(/r) is less than one. In order to ensure this, we perform a 2-stage Greedy 
similar to e.g. [7, 10, 13, 14, 15]. In the first step, we determine a preliminary snap¬ 
shot set Sm such that < 1 for ah € Strain- This is realized in lines 2-7 in 

Algorithm 2. The second loop in lines 9-14 enriches the so determined Sm so that for 
the resulting set S'at, we get max^gStrain "A ^toi- In addition, we determine an 

orthonomal set of functions and define the RB trial space as A^r := span(AAr). 


Algorithm 2 Greedy Algorithm 

1: choose G X> arbitrarily, ■<— m*’^(^i), A ^ 1, Sn '■= Xn '■= {Ci} 

2 : while max^gHtrain 'rAf(Ai) ^ 1 do 
3: ^lN+l ^ argmax^g2^_,^._^ Tjv(/r) 

4: Sn+1 S'tv U {flN+l}, iN+1 <— (flN+l) 

5: orthonormalize Cn+i w.r.t. Xn — ^n+i, Xn+i t— Xn U 

6: N-^N + l 

7: end while 
8: M 

9: while max^gHt„i„ An(p) > £toi do 
10: ^N+1 <- argmax^gStrain 

11: Sn+ 1 t— Sn U {^jv+i}i iN+i t— u*’’^{^ n+i) 

12: orthonormalize ^at+i w.r.t. Xn —t Cat+I: ^n+i t— Xn U {^at+i} 

13: N ^ N + 1 

14: end while 


5. Numerical Experiments. We now present results of some numerical exper¬ 
iments. To this end, we use the following data: T = 1, 12 := (—150,150) C M, d = 1, 
V = [0,100]. In view of Remark 2.1, we consider quadratic abatement costs 7 (t, x)"^/2 
which are discounted to the end T of the trading period with an interest rate of 0.05. 
This results in the following running costs: 

(r(M))(t,x) 

As a (single) parameter /i G R, we chose the overall amount of emissions by 
the involved companies that would arise without any reduction motivated by the EU- 
ETS, i.e., the “business as usual emissions”. The control 7 (t, x) denotes the amount of 
avoided emissions so that ^ — x) is the remaining emission, i.e., the drift. Together 
with a constant diffusion term, this yields: 

(A^(/x; u)) := x) - x) - x) 
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The permits expire worthless at the end of the trading period and the penalty pay¬ 
ment is normalized to one currency unit, which is reflected by the terminal condi¬ 
tion ut{x) := x~^ for a; S fi in (2.8c). The Neuman conditions are modeled by 
gieit{t) = 5ieft = 0 at a; = —150 and 5right(i) = ffright = 1 at x = 150. The function 
7 M- L^{g,;u) := dtu + A'^{ijl\u) is so that the optimal control can in fact be 
characterized as a critical point. 

5.1. Discretization. The discretization is done by finite differences, fully im¬ 
plicit in time with step size At = and using central differences on a regular grid 

with mesh size Ax = 1.5 = In such a discretization, the discrete representa¬ 
tion of the derivatives is a linear combination of the value function u at different 
discretization points. 

The obvious fact that the above defined u) and /"’'(p-) are smooth in 7 and 
u thus also holds for their discretizations. Due to this smoothness, the assumptions 
of Lemma 3.3 hold. In order to stabilize the discrete equations, we use artificial 
diffusion which is chosen parameter-independent in order to obtain a /i-independent 
discretization. Clearly, A"*{fj,; u) and also are separable in the parameter which 

also translates to their discrete versions. 

5.2. Estimation of the inf-sup constant. We start by reporting results con¬ 
cerning the computation of the lower bound in (4.5) for the inf-sup constant 

using the anchor point strategy described in Algorithm 1 above. We compute an 
‘exact’ value of by solving the corresponding high-dimensional (‘truth’) gen¬ 

eralized eigenvalue problem. The results are shown in Figure 5.1. The thin lines 
represent the intermediate lower bounds which are generated by successive iterations 
of Algorithm 1. Different colors represent different stages of the algorithm. The final 
lower bound on the training set is marked with the thick line. 


0.2 

0 

- 0.2 


-0.4 


- 0.6 


- pN{g) 


0 20 40 60 80 100 

iV = 10 M 



Fig. 5.1: Lower bound PYP{g) for the inf-sup constant Pnig)- Different colors indicate 
iterations of the algorithm. 


5.3. Error estimator. Next, we test the performance of the error bound Ajv(/i) 
and the indicator r)0®(^) given by the BRR theory. 

To obtain the Lipschitz constant p, we use the decomposition given by Lemma 
3.3. Since / depends quadratically on 7 and L depends linearly on 7 and on v for the 
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EU-ETS, it holds that Q 2 = 92 ~ 0. Furthermore, 


= sup 

7^ gRiiAT 


£>0 = sup 

gRdAT 


57 /- 


A/',7^ 




IIt^IIr-^^ 

||L(RAf,RAf) 


= sup 

7-^ gRdAT 


h' 


Xp0.05(t-T)| 


L(R'i,R)J^ 


h^ll 


117^11 


RdAf 


, and 9 ^ = sup 

a;Af 


RdAf 

yL^U^l 


L(R‘i,R)J^ 


l|a;^llx^ 


which can be computed solving generalized eigenvalue problems. Since these quanti¬ 
ties are independent of /r, is suffices to determine them once in the offline phase. 

The norm of the residual ||^^(/r)(x^(/r))||YV can be computed efficiently in the 
online phase using the separation in the parameter. Inserting (4.2) and xn{p-) = 
Yl,n=l {9.)^n into (/x)))yv results in a sum of products 

of /i-independent functions which can be evaluated fast in the online phase and Y-^- 
dot-products which are precomputed in the offline phase. 

In Figure 5.2, we compare the true error (again computed with respect to the 
detailed discretization) with the indicator, the error bound and the size of the residual. 
The values correspond to the maximum of the corresponding quantities over a test 
sample in V. We recall that the BRR-based error estimator is only meaningful for 
T < 1, which is here the case for N > 10. 



-T)Y®(/i), indicator 

-A 7 v(lt), error bound 

-true error 


-residual 


Fig. 5.2: Error, indicator, bound and residual over N. 


In order to investigate the effectivity of the error bound, we fix A = 10 and 
vary the parameter 9 over the parameter range 27 = [0,100]. We can see in Figure 
5.3a that indicator, bound and residual are numerically zero for the snapshots, where 
the true error of course also vanishes. The effectivity of the error bound, i.e., the 
ratio of error estimator and true error, is shown in Figure 5.3b. We expect a growth 
with respect to increasing 9 , since the PDE starts becoming increasingly convection- 
dominated. However, the maximum size is below 8 which seems a reasonable size to 
us. Of course, well-known techniques for stabilization as well as parameter-adaptivity 
may additionally be used, e.g. [9, 11]. 

6 . Summary. We have presented a Reduced Basis Method (RBM) for rapidly 
solving the parameterized Hamilton-Jacobi-Bellman (HJB) equation with the specific 
application to the European Union Emission Trading Scheme (EU-ETS). In particular, 


15 









(a) Error, indicator, bound and residual over (b) Effectivity of the error bound. 

^ T). 


Fig. 5.3: Effectivities over parameter range T> = [0,100] for N = 10. 


we have introduced a rigorous bound of the error with respect to an online-efficient er¬ 
ror estimator. The involved parameter-dependent constants can be computed online- 
efficient by an anchor-point based Successive Constraint Method. Numerical experi¬ 
ments confirm the effectivity of the estimator. 

Future research will be devoted to the specific application of the HJB-RBM to the 
EU-ETS in order to determine optimal regulatory strategies. From the mathematical 
point of view, we will consider extensions to more general settings, i.e., more general 
nonlinearities and stronger convection in the operator. 
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